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We present a worm sampling method for calculating one- and two-particle Green’s functions us¬ 
ing continuous-time quantum Monte Carlo simulations in the hybridization expansion (CT-HYB). 
Instead of measuring Green’s functions by removing hybridization lines from partition function con¬ 
figurations, as in conventional CT-HYB, the worm algorithm directly samples the Green’s function. 
We show that worm sampling is necessary to obtain general two-particle Green’s functions which 
are not of density-density type and that it improves the sampling efficiency when approaching the 
atomic limit. Such two-particle Green’s functions are needed to compute off-diagonal elements 
of susceptibilities and occur in diagrammatic extensions of the dynamical mean field theory and 
efficient estimators for the single-particle self-energy. 

PACS numbers: 71.27.+a, 02.70.Ss 


I. INTRODUCTION 

The Anderson impurity model (AIM) 1 ^ is one of the 
fundamental models for electronic correlations. The 
model was originally developed to describe the physics of 
magnetic impurities in solids, but nowadays also serves as 
a model for quantum dotspEl adatoms on surface^ 7 and 
appears as an auxiliary model in the context of dynamical 
mean field theory (DMFT)PED Continuous-time quan¬ 
tum Monte Carlo (CT-QMC) algorithms'^ are state 
of the art for the numerical solution of the AIM. They 
are based on a stochastic sa mplin g of an imaginary time 
partition function expansion. 17 18 

The methods are formally numerically exact and, in 
contrast to other impurity solvers CsHSJ?] ca n treat impuri¬ 
ties with many degrees of freedom, general interactions, 
and continuous bath dispersions. The most widely known 
representatives are formulated as an expansion of the par¬ 
tition function eithe r in terms of the interaction (CT-INT 
and CT-AUX) 12 ! 15 ' or in ter ms of the impurity-bath hy¬ 
bridization (CT-HYB) JISMi with the resulting series sam¬ 
pled stochastically. 

A variant of continuous-time algorithms, usually re¬ 
ferred to as the worm algorithm, expands both the par¬ 
tition function and the Green’s function. This results 
in the configuration space sampled by Monte Carlo to 
be enlarged (see Fig. [l]). This concept has been pio¬ 
neered for diagra mmat ic Monte Carlo solvers for bosonic 
Green’s function d 18 ) 26 ! and adapted for fermionic one- 
particle Green’s functions for the CT-INT algorithm. 27 

In this paper, we introduce a generalization of the 
worm algorithm for the (multi-orbital) hybridization 
expansion. 14 While worm sampling is not restricted 
to any specific quantity, we show the application to 
fermionic two-particle Green’s functions which are neces¬ 
sary to compute response functions and which appear in 
formulations of non-local extensions of the DMFT such as 
the dynamical vertex approximation^^) the dual fermion 



FIG. 1: Illustrating the concept of worm sampling. The con¬ 
figuration space of the partition function Cz is enlarged by the 
configuration space of the n-particle Green’s function C G („). 
A random walk in the combined configuration space is shown, 
where dashed lines represent the transition moves between the 
two configuration spaces and solid lines the moves within one 
space. 


approach,^ the one-particle irreducible approach 30 and 
the DMFT to functional renormalization groupPb) They 
also appear in the measurement of single-particle self¬ 
energies using the ‘improved estimator^ 2 technique that 
has been shown to yield high precision estimates for the 
high-frequency behavior of Green’s functions. 

In Section [TT] we motivate of our work by showing that 
conventional CT-HYB partition function sampling fails 
due to ergodicity problems when approaching the atomic 
limit and when calculating general two-particle Green’s 
functions. 

Section |III| first gives a short overview of worm sam¬ 
pling and then generalizes CT-HYB to the Green’s func¬ 
tion space, introducing the Monte Carlo update proce¬ 
dure of our CT-HYB worm method. Section IIVI intro¬ 
duces the measurement procedure. Section [V] presents 
the results for large interactions and the atomic limit, 
where analytical solutions are available. Section VI fo¬ 
cuses on results for the two-particle Green’s function of 
the two-orbital model, further validating the worm sam¬ 
pling algorithm by exploiting the SU(2) symmetry of the 
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magnetic (spin) susceptibility, 
brief summary. 


Section VII provides a 


II. MOTIVATION 


We start with a brief motivation for measuring the 71- 
particle Green’s functions G^ with worm sampling. The 
Hamiltonian considered here is that of the multi-orbital 
AIM, which in its most general form reads: 


-Haim — ^ ^ ( Uap^sdil^d^dsd^ V ^ ) - a d) v d ri T 
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Here, 4 {d a ) denotes the creation (annihilation) opera¬ 
tor of an electron with spin-orbit flavor a on the impurity 
and 4 0 (c fca ) denotes the creation (annihilation) opera¬ 
tor of an electron of momentum k in the non-interacting 
bath that belongs to the impurity flavor a. The impu¬ 
rity problem is characterized by the one-particle levels 
i a , the interaction matrix U a p^s, the bath dispersion 
Ska, and the hybridization strengths . In CT-HYB, 
given some inverse temperature /3, the partition func¬ 
tion Z = Tr e~P HA1M of the AIM Hamiltonian (|T]) is ex¬ 
panded in terms the hybridization Hhyb- The trace then 
decouples into a local part described by the local impu¬ 
rity Hamiltonian (H\ oc ) and a bath part described by the 
conduction electron bath (Hbath)- With the bath parti¬ 
tion function Zbath = Tre~ /3Hhath we find (for a detailed 
derivation see Ref. [13]): 



dr k -i 

-2 



x Tr 


T T e /3Hloc da fc (T fc )4 fc _ 1 (' r fc- 1 ) - da^dl^n) 


= W loc (k,T 1 ,...,T k ) 


x detA, (2) 

=Wbath(fc,Tl,...,Tfc) 


Here 4(4 ( d a (r )) are the operators of Eq. (JlJ in 
Heisenberg representation, whose evolution in imaginary 
time is given by H\ oc . Further, T t is the Wick’s time or¬ 
dering operator and A denotes the ( k/2 ) x {k/2) matrix 
of all possible hybridization lines between T\...r k , with 
the elements: 


k, 7 


(vrrK' 1 

e @Ek-y _j_ \ 


— e -£k-y(T-0) 


T > 0 
T < 0 


(3) 


where r = Tj — Tj. We refer to hybridizations as diag¬ 
onal if A aa'(r) = 0 Ma ^ a'. Otherwise we call them 



FIG. 2: Illustration of a configuration of the partition function 
Z for an expansion order of k/2 = 3. Here we show the case of 
a flavor-diagonal hybridization function, which only connects 
operators of the same flavor to one-another. The different 
flavors are denoted using different colors (red, blue). When 
connecting creation (filled shapes) and annihilation (empty 
shapes) operators by hybridization lines, we number all cre¬ 
ation operators from 1 to k/2 and all annihilation operators 
from 1 to k/2. 


off-diagonal. In the following we will restrict ourself to 
diagonal hybridizations, even though in principle also off- 
diagonal hybridizations can be considered. 

In Eq. § abbreviations for the local weight wi oc and 
the bath weight Wbath are introduced, which become im¬ 
portant when defining the Monte Carlo algorithm. In 
Fig. [2] we provide an illustration of a configuration of the 
partition function Z for a given expansion order k/2. A 
more detailed discussion is found in Ref. m- In this 
work we refer to the expansion order as k/2 so that the 
number of operators in the local trace is given by k (dif¬ 
ferent conventions exist in literature). 

When measuring the one-particle Green’s function 
G«(t) in conventional CT-HYB, sampling takes place 
in partition function space Cz, i.e., the orders k/2 and Tj 
in Eq. are sampled. One starts from the functional 
identity: 


^ Ctot.' 


(4 


1 sz 

Z 5A a > a (—T)' 


(4) 


The conventional estimator is obtained from Eq. ([2]) by 
replacing the functional derivative in Eq. Q with the 
partial derivative and using the chain rule, which gen¬ 
erates local operators by detaching their hybridization 
lines (a detailed derivation can be found in Appendix C 

of Ref. 1551 ) : 


G&(r) 


1 

0 



nm 


det A (ram ) 
det A 


x sgn • S(t, T m - T n )5 aam S a 'a n ) , (5) 

/ MC 


where A^™") is the {k/2) x {k/2) hybridization matrix A 
with the n-th row and m-th column removed, correspond¬ 
ing to the removal of hybridization lines; S{t, r m — T n ) 
specifies the imaginary-time bin of the measurement and 
(...)mc refers to the Monte Carlo expectation value of the 
r integrals and k sums of Eq. Q including the weighting 
factor e ~P Hloc \ “sgn” denotes the sign imposed by the 
Wick time ordering. In the following we will denote the 
estimate Eq. M, suppressing the indices a , a', as G^(r). 
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Computing the Green’s function by evaluating the quo¬ 
tient of the hybridization matrix reveals a first short¬ 
coming of this approach: the estimator in Eq. © fails 
if the hybridization between the impurity and the bath 
becomes very weak. We will later see that the estimator 
of worm sampling instead does not depend on the deter¬ 
minant ratio of the hybridization matrix A, such that 
sampling is still possible for small or vanishing A. This 
suggests that worm sampling does a better job for sys¬ 
tems approaching the atomic limit. We point out that 
some methods exist in order to improve the estimator 
in Eq. ([5]). A recent approach is the so called remove- 
shift measurement (or sliding measurement), which has 
been implemented for density-density codes.®! While the 
remove-shift estimator is capable of enhancing measure¬ 
ments by decreasing auto-correlation times, it does still 
depend on operator pairs which are connected to the bath 
over their hybridization. As such, this approach does not 
cure the problem encountered for weakly hybridizing sys¬ 
tems. 

Eq. @ is restricted to diagrams produced by parti¬ 
tion function sampling and does not generate off-diagonal 
Green’s function contributions for diagonal hybridization 
matrices A aa . However, for non-density-density interac¬ 
tions, such terms are indeed present in the two-particle 
Green’s function G^ 2 K One can immediately see this 
for the SO(n) SU(2)-conserving Slater-Kanamori in¬ 
teraction: here, the spin susceptibility is invariant un¬ 
der spatial rotations, such that, e. g., (5' i: (r)S' 2 (0)) = 
(5'a:('r)<5 a ;(0)). The spin susceptibility in ^-direction re¬ 
lates to flavor-diagonal terms of G^: 

{Sl(r)Si{ 0)) = |((r4(r) - nl(r))(n J t {0) - n^(0))> = 

^( c lt( r ) c it(' r ) c jt(°)Gt(°) - 4r( T ) c q-( T ) c ii(°)ot(o)- 

c li( r )^t(T)S't(°) c tt(0) + c^(r)c^(r)c^(0)c^(0)}. (6) 

All terms can be obtained in conventional CT-HYB by 
removing one hybridization line for orbital i and one for 
orbital j, analogous to Eq. The spin susceptibility 
in x-direction on the other hand manifests itself as spin 
flip terms in G^ 2 \ which are off-diagonal: 


{SUr)Slm = i((4(r) + SUr))(Si( 0) + SL( 0))> = 

|<4t( T ) c it( T )4r(°)Gt(°) + c lt( r ) c it( T ) c ]4.(°) c Jt(°)+ 

cU r ) c *tMc}t(°) c tt( 0 ) + c^(r)c it (r)c^(0)cj 1 .(0)>. (7) 


We emphasize that the two-particle generalization of 
Eq. ([5]) does not provide the spin-flip terms of Eq. ([T]) 
but only density-density-like terms as in Eq. ©>■ One can 
obtain Eq. ([ 7 ]) by a functional derivative as in Eq. (4j), 
albeit with a hybridization function A aa i that is either 
off-diagonal in the orbitals or the spins. Such terms are 
however not generated in the hybridization expansion 
Eq. @, at least not for an orbital- and spin-diagonal 
A™™. 


A particularly important application of these off- 
diagonal elements is found when extracting the self¬ 


energy E(zw) from the equation of motion; a technique 
that leads to precise high-frequency estimates. This 
method is usually referred to as improved estimators and 
has so far only been implemented for density-density 
interactions.^ For interactions of non-density-density 
type, off-diagonal terms of the two-particle Green’s func¬ 
tion are needed when implementing improved estimators 
for the self-energy and the reducible vertex. We will 
show how worm sampling is capable of supplying such off- 
diagonal terms, hereby overcoming the systematic short¬ 
coming of traditional CT-HYB algorithms of being re¬ 
stricted to Green’s functions generated by the type of 
AIM hybridization. 


III. SAMPLING AND ERGODICITY IN 
GREEN’S FUNCTION SPACE 


In order to solve the restrictions of the conventional 
Green’s function estimator, Eq. <©, we may be tempted 
to turn the diagrammatic series of a local observable O 


(0(r)) = 


-'bath 


E 


I dr k / dr k -i ... / dr\ x 
ke 2N 0 ,a fc •' Tfe — 1 -' r f=-2 Jo 

x TY[T T e~ ljHl ° c O(T)d ak (T k ) ■■■ d a2 (T2)d} ai {Ti)] 

x det A (8) 

into a Monte Carlo estimator by inserting O into dia¬ 
grams from the expansion of Z 1 Eq. ([ 2 ]), and measuring 
the weight ratio. However, as already noted in Ref. [161, 
such an estimator for the Green’s function is not ergodic 
(we will elaborate on this in Section III BI. 

Using worm sampling, we solve this issue by enlarging 
our configuration space 


C = Cz © C G (to 


(9) 


to contain both types of diagrams of Eq. ([2| of the par¬ 
tition function space Cz and Eq. ([8]) of the ?r-particle 
Green’s function space C G („) (see Fig. [l]). The sampling 
in C G (t>) allows us to generate all diagrams for the Green’s 
function, thereby circumventing the ergodicity problems 
of both the estimator constructed from insertion of local 
operators and from removal of hybridization lines. While 
C G (n) was originally introduced as an auxiliary space to 
restore ergodicity and lower auto-correlation times for 
Cz, here the reverse can be argued: excursions to parti¬ 
tion function space lower the auto-correlation times and 
provide the proper normalization for the Green’s function 
(cf. Section IV). 

In this work we restrict ourselves to sampling as 0(t) 
in Eq. (JsJ) the one-particle Green’s function and the two- 
particle Green’s function in imaginary time r defined by: 


GaE (n, tj) = —(T T d ai (n)dl 2 ( tj )) (10) 

ct±oi20i3CX4 (A , Tj , T k , T;) 

{T r d ai (T i )dl 2 (T j )d a3 (T k )dl ti (T l )} (11) 
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(a) 

(c) 


P\ - (j)- 
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over all degrees of freedom of the Green’s function G 


= J dr l- dr « G i" ) 1 ...,an( T l)- 7 V*)- ( 12 ) 

CKl ,...,cx. n 
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(f) 
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(e) 
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This is not a “physical” partition function in the sense 
that it is connected to a thermodynamic potential, but it 
simply represents a phase space volume in Green’s func¬ 
tion space. We will now discuss all the steps mentioned 
in Fig. [3] in full detail. 


A. Worm Insertion and Removal Steps 


FIG. 3: Illustration of Monte Carlo moves in the extended 
configuration space of worm sampling. Circles denote oper¬ 
ators connected by hybridization lines (indicated by vertical 
lines), while rectangles denote worm operators. Moves (a), 
(b) and (c) correspond to insertion, removal and shift of an 
hybridization operator pair in the configuration space C z , re¬ 
spectively. Labels (d) and (e) exemplify worm insertion and 
removal moves transitioning between the two spaces, where 
the parameter j/b rescal es the phase space volume of C G (i) 
[for details see Sec. |III A . Labels (f) and (g) denote removal 


and insertion of an hybridization operator pair in C G ( 1 ) [Sec. 
|III B| ; (h) labels the worm operator replacement move in C G ( 1 ) 
|Sec. IlIICl. 


Restricting worm sampling to the Green’s functions 
space Cq(tl) has two reasons: (i) the one and two-particle 
Green’s functions include almost all relevant information 
about the quantum impurity (see Eq. § -Eq- ©)■ (h) 
when sampling the one- and two-particle Green’s func¬ 
tion, we can compare our results against the measure¬ 
ments in the partition function space Cz (especially with 
regards to the normalization, error bars and strong in¬ 
sulating cases). While a similar comparison in principle 
would be possible for the three-particle Green’s function 
G< 3 \ we do not consider this quantity because of the 
high computational effort involved and the less physical 
significance in comparison to G^ and G^ 2 '. 

In Fig. [3] the Monte Carlo moves in Cz and C G („) are 
illustrated. We included all steps needed to be ergodic 
and to decrease auto-correlation lengths in both con¬ 
figuration steps. The pair insertion and removal steps 
in Cz (Fig. [3|a),(b)) are typical in the CT-HYB algo¬ 
rithm. We further introduce the operator shift move for 
Cz (Fig. [3|c)), which shifts the time of a creation or an¬ 
nihilation operator. 

For later discussion, we set up a modified partition 
function Z G ( n ) in configuration space C G (n) by integrating 


The worm insertion and removal steps are transition 
steps between the two configuration spaces, depicted in 
Fig. [ 3 ] (d),(e). In order to sample in Cz and C G <», jump¬ 
ing between the two spaces is needed. In general, the con¬ 
figuration spaces Cz and C G ( n ) have very different phase 
space volumes. This difference is balanced out by intro¬ 
ducing a weighting factor 7 /") so that the total partition 
function reads 


W = Z + »?<"> Z G( »). (13) 

For now it was not formalized how rf n ^ scales with 
the number of orbitals, temperature and interaction 
strength. It is best to choose rf n ' 1 so that the simula¬ 
tion spends an equal amount of steps in Cz and C G ( n ) . 
We revisit this fact when discussing the normalization of 
the worm result in the following section. 

It is important to mention that the only difference be¬ 
tween worm operators and hybridization operators is the 
missing of hybridization lines. This has some implica¬ 
tions for our Metropolis acceptance rates. The proposal 
rate of inserting a worm is given by the same expression 
as the proposal rate of inserting n hybridization operator 
pairs, i.e.pS 


f(Cz^C GM )= d ^. (14) 

Adding worm pairs results in the expansion order k/2 
of the local trace being increased by n, whereas the ex¬ 
pansion order in the determinant is kept constant. This 
adds an ambiguity to the expansion order which needs to 
be kept in mind. The weight of a configuration in 
modified by r/") is then: 

P(C G (.n ), Ti , ... , TfcJ Tjj , ..., Ti 2n ) = 

?7 (n) • wioc (k + 2n,n,... ,Tfc;Ti 1} ...,Tj 2 „) 

tt>bath(fc,Ti, ... ,T k )dTi ... d,Tk- (15) 

We point out that combining the proposal probability 
and the configuration of the weight, the 2 n infinitesimals 
drq ... d,Ti 2n do not cancel as they would have in partition 


























5 


function sampling. This is due to the extra local degrees 
of freedom introduced by the worm and is integrated over 
in the computation of Z G („) (12). The proposal probabil¬ 
ity for removing the worm is simply: 


f{C G (n) —> Cz ) — 1- (16) 

Note, since there is only one worm in the trace at a given 
time, we always propose to remove exactly this worm. 
The Metropolis acceptance rate of a worm insertion is 
hence: 



a(Cz C G (n )) — 

' („) |wi oc (fc + 2n, Ti, ... ,r k > Ai i i Ft 2Tl )| o2n 

|wi oc (fc,Ti, ... , Tfe)| 


(17) 


The bath weight iCbath, which includes the hybridization 
matrix, cancels out due to the fact that the bath remains 
unchanged. 

The inverse gives the acceptance probability of a worm 
removal: 


FIG. 4: An “insertion estimator”, i.e. the mere insertion 
of local operators into a diagram from Cz without sampling, 
is not ergodic: it fails to produce diagram (3) because (2b) 
violates the Pauli principle and is therefore never reached. By 
first transitioning to C G ( i) space from (1) and then inserting 
a hybridization operator pair into (2a), one indeed is able to 
reach diagram (3). 


Acceptance rates are similar to the corresponding accep¬ 
tance rates in Cz space: 


a (C G ( n ) —> Cz) — 


min 1, 


JL_ |wioc(fc,T~ 1; ,T k )\ _ 

Hoc(fc + 2n, Ti,... ,T fc ; r il ,...,r» 2 „)| 


1 

P2n 

(18) 


a(CQ( n ) 5 k T 2 n —^ k T 2 n -{- 2) — 

|wioc (k + 2 n + 2 ,Tl, ... jTij,... ,Ti 2 „, ... , Tk',Ti,Tj) 


|«'ioc(fe + 2n, n 5 ••• ? 7~ii j ••• j T~i 
|Wbath(fc + 2, Ti, ... , Tk', Ti,Tj)\ /3 

\Wba,th{k,Tl, ... ,Tk)\ {{k + 2)/2) 2 


1'2n ? ”• > 
2 


(19) 


We point out that we jump between Cz and C G o> and 
between Cz and C G ( 2 ), but never between C G ( 1 ) and C G < 2 ). 
As mentioned in Section|TTJ the two-particle Green’s func¬ 
tion for non-density-density interaction includes spin flip 
and pair hopping terms. The one-particle Green’s func¬ 
tion, on the other hand, is always flavor-diagonal for 
flavor-diagonal hybridization functions. This way, insert¬ 
ing two worm pairs consecutively by attempting to jump 
from Cz to C G ( 1 ) and then to C G ( 2 ) will fail to provide 
flavor-off-diagonal i.e. spin flip and pair hopping terms. 
A very similar observation was recently made for the con¬ 
ventional CT-HYB algorithm with a flavor-off-diagonal 
hybridization function 35 . 


where the worm operators are located at times 
Tjj,... ,Tj 2n . The Metropolis acceptance rate for a pair 
removal in the Green’s function space is then just given 
by the inverse of Eq. (19). 

We remind the reader of the fact that the local weight 
w\ oc in Eq. (19) is expressed relative to a factor k+2n+2 , 
while the bath weight u>bath is expressed relative to a 
factor k + 2. The discrepancy comes from the n worm 
operator pairs in the local trace without hybridization 
lines. 


C. Worm Replacement Step in Green’s Function 
Space 


B. Pair Insertion and Removal Steps in Green’s 
Function Space 

In order to generate all possible Green’s function con¬ 
figurations, we need to introduce additional updates in 
the Green’s function space C g m . This is a crucial part of 
worm sampling: without it, the estimator is not ergodic 
(cf. Fig. [4]). 

This explains why we are required to sample the 
Green’s function space C G { n ) separately with operators 
having hybridization lines attached. To this effect, we 
perform insertions and removals of hybridization opera¬ 
tor pairs also in Green’s function space (Fig. |3](f),(g)). 


While insertion and removal moves formally fulfill 
the condition of ergodicity, worm sampling requires a 
shift/replacement move in order to allow for acceptable 
auto-correlation lengths. We elaborate on this require¬ 
ment here. 

Let us assume a local trace filled with hybridization 
operator pairs. We now attempt to insert a worm pair 
into this trace. It turns out that inserting a worm pair, 
where the worm operators are relatively close to one an¬ 
other is probable, while inserting a worm pair where the 
worm operators are far apart is less probable. This is 
because of (i) possible quantum number violations since 
there may be many creation and annihilation in between 
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the pair for long time differences, and (ii) the pair inser¬ 
tion might lead to an energetically disadvantageous local 
configuration which is unfavorable to have for a long time. 

Problem (i) is especially severe if we have a large 
amount of operators in the trace, which occurs at small 
interaction or low temperatures. Additionally, more re¬ 
strictive interaction types, such as the density-density in¬ 
teraction, produce more rejects due to quantum number 
violations of attempted worm inserts. This is why we do 
not observe this auto-correlation problem at high tem¬ 
peratures, high interaction parameters and more general 
interactions such as Slater-Kanamori interactions (which 
may change the quantum number in the local trace). 

The solution to this problem is found in 
shift/replacement moves. We consider, instead of a 
general worm shift move, a replacement move which 
exchanges one of the worm operators with an operator 
of the hybridization expansion, i.e., we replace it with 
one of the same flavor connected by a hybridization line 
as illustrated in Fig. §h)- 

This way we do not have to recalculate the local trace, 
as two locally indistinguishable operators switch position. 
Instead, we need to recalculate the determinant of the hy¬ 
bridization matrix since the replacement corresponds to 
a shift of the worm operator and a shift of the hybridiza¬ 
tion operator. Further we do not encounter any rejects of 
proposed moves due to local quantum number violations. 

It turns out that worm replacement moves (or in the 
same way worm shift moves) are equally important for 
traces with very few operators because of problem (ii). 
This problem typically occurs if the weight e~ UT of the 
worm becomes prohibitively small, i.e., in particular for a 
large interaction strength and a long r difference such as 
|. We are then effectively restricted to inserting operator 
pairs into the trace, which are very close to each other in 
imaginary time. These pairs have similar properties as 
density operators and can in principle be inserted for very 
high insulating cases. By inserting hybridization pairs at 
short distances T t — Tj and then replacing one worm oper¬ 
ator with one hybridization operator we are able to pass 
this restrictions of the time evolution. As we will show 
in the following, the replacement move only depends on 
the ratio of the determinant of the hybridization matrix. 

The proposal probability of a worm replacement step 
is given by: 

f'(C G (n ), k + 2n —> k + 2n) = . (20) 


This corresponds to selecting one creation/annihilation 
operator of the 2 n worm operators at random and select¬ 
ing one creation/annihilation of the same spin-orbit fla¬ 
vor with a hybridization line. In practice, we choose an 
operator from the k/2 operators of the same type (anni- 
hilator/creator) and then discard flavors, which are not 
equivalent to the worm flavor. The proposal probability 
of switching the operators back to their original position 
is hence also given by Eq. (201. 



0 20 40 60 80 100 120 140 160 180 200 


T 

FIG. 5: One-particle Green’s function G^\t) in imaginary 
time r, illustrating the ergodicity problem of the worm algo¬ 
rithm for an average expansion order of k/2 ~ 40. Param¬ 
eters: inverse temperature /3 = 200 /D, Coulomb repulsion 
U = 0.5 D and n = 0.3 D (out of half-filling) for the the single¬ 
orbital AIM with semi-elliptic conduction electron density of 
states with half-bandwidth D = 1 and V = 0.5 D. The bal¬ 
ancing parameter r/ 1 ) was chosen in the interval [0.15,0.22]. 
We observe the ergodicity problem between r = 25 /D and 
r = 175 /D (blue curve). When adding replacement moves, 
we are able to insert worm operators for such r’s around f3/2 
(green triangles) and hence obtain much better results. We 
have additionally supplied G (1 > (t) for the measurement in 
partition function space (red curve). 


We observe that the proposal probabilities for the re¬ 
placement move cancel out and the acceptance ratio is 
fully determined by the ratio of weights. Further, the 
local weights cancel, since a worm operator and the cor¬ 
responding hybridization operator are indistinguishable 
within the local trace. The Metropolis acceptance rate is 
hence given by: 

a (C G ( n ), k T 2 n — 1 k T 2 ti ) = 

“A 1 ' |«W*,n.*,. <21) 

where t* refers to the initial position of the worm opera¬ 
tor and Tj to the initial position of the operator with the 
hybridization line. Fig. [5] shows how worm replacement 
moves alleviate the ergodicity problem of the worm algo¬ 
rithm for the situation where many operators are found 
in the local trace. 

We would like to use the opportunity to point out the 
difference between a worm replacement and a worm shift 
move. The acceptance rate of the worm replacement 
move depends on a determinant ratio of two matrices 
of dimension (k/2) x (k/2), where k here refers to the 
number of operators with hybridization lines connected. 
In that sense it is very comparable to the determinant 
ratio of two matrices of dimension (k/2 — 1) x (k/2 — 1) 
and (k/2) x (k/2) in Eq. § when changing the order be¬ 
tween k/2 and (k/2 — 1) in partition function sampling. 
The acceptance rate of a worm shift move, on the other 
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hand, only depends on the ratio of the local traces. While 
for the worm replacement move we are able to pass the 
restrictions of the local time evolution, for the worm shift 
move we are able to pass the restrictions of the hybridiza¬ 
tion function. When calculating strong insulating cases 
we profit the most if we consider both moves. 


IV. WORM MEASUREMENT 


We now show how the measurement of Green’s func¬ 
tion looks in Cqm . It turns out that the measurement 
itself is trivial and we only need to find the correct nor¬ 
malization of the Green’s functions measured and the 
correct sign. For the one-particle Green’s function G ^ 
a worm is defined by the operators d(rj) and <v(rj). The 
correct weight is intrinsically given as we sample in the 
Green’s function space C G („). Thus, the estimator of the 
Green’s function simply follows as: 

42 ( t ) = ( s § n ' 5 ( r > T i ~ t j))mc- (22) 

The Green’s function in Matsubara frequencies can be 
calculated by substituting the d-function by the Fourier 
transform: 

GgM = ( sgn-e^-^MO (23) 

The measurement of the two-particle Green’s function 
in Matsubara frequencies in the particle-hole channel is 
given by: 


42 (* v i * w ) = 


(sgn • e i ^-Tj) e ^ , (r k -T l ) e Mn- T i)) Mc _ 


(24) 


The imaginary time arguments 7 $,..., 77 are assigned to 
creation and annihilation operators according to Eq. (113- 
While we both employ Eq. (22) and Eq. (231 for the 


one-particle Green’s function measurement, the measure¬ 
ment of the two-particle Green’s function in Matsub¬ 
ara frequencies, Eq. (24), is far more convenient than 


a binned measurement in imaginary time. It is especially 
difficult to resolve jumps in the imaginary-time measure¬ 
ment due to fermionic sign changes in the time ordering of 
operators. Measuring the two-particle Green’s function 
in imaginary time using a binning procedure and then ap¬ 
plying the Fourier transform gives wrong high frequency 
asymptotics, while the direct measurement in Matsubara 
frequencies is free of errors resulting from binning. 

As with conventional sampling, we do not observe any 
sign-problem for worm sampling in the case of a flavor- 
diagonal hybridization function. However, unlike in the 
G C J estimator, the flavor indices and the imaginary time 

bins in the worm estimator 42 are outer indices, such 
that the mean sign in principle also becomes flavor and 
r dependent. 


Eq. (22) and Eq. (231 are normalized to Z G (e >, Eq. (24) 
to Z G ( 2 ) as defined in Eq. (12), as opposed to the phys¬ 
ically correct normalization to Z. We will now discuss 
the normalization in more detail. 


A. Normalization and Auto-Correlation 


In principle we are ergodic in C G <„) , when assuming 
worm replacement or worm shift moves. It turns out 
however that we need to sample both in C G( „) and Cz 
with about the same number of steps to fix the normal¬ 
ization z °f the thermal expectation value in Eq. |8j). 

When measuring the Green’s functions in C G („> we 
implicitly normalize with the number of steps taken in 
C g m- We correct for this factor by explicitly counting 
how many steps N G were taken in C G (n) . We further 
count how many steps Nz were taken in Cz- This esti¬ 
mates the size of the configuration space Cz , which then 
gives the correct normalization. The normalization for 
is then given byp3 


Q ( n ) 


1 N g („) 

7y(n) Nz Ca ’ 


(25) 


where G^J is measured in C G ( and the factor l/rj^ is 
a result of rescaling Z G ( n ) in Eq. (13). 

Let us note that Eq. (251 is only one way of normaliz¬ 
ing the worm measurement. In a different approach, we 
could do the entire sampling in worm space, without re¬ 
moving the worm operators at all. We are then required 
to generate worm configurations by shift moves and re¬ 
placement moves. In this case, we could normalize the re¬ 
sult by assuming some physical knowledge of the Green’s 
function. One possibility is to extract the normalization 
by assuming the correct behavior of the large-frequency 
asymptotics of G G \iu) or G^ (iu, iv', iui). 

In order to calculate the Monte Carlo expectation value 
Eq. (23), we still need to divide by the number of mea¬ 
surements N taken. It is important to notice the differ¬ 
ence between the number of measurements N and the 
number of steps N G and Nz taken since it is common 
to skip steps during two consecutive measurements to 
assure uncorrelated measurements. 

This directly relates to the auto-correlation length of 
the QMC sampling. The auto-correlation length in worm 
space C G (n) looks very different from the auto-correlation 
in partition function space Cz- A well-accepted estimate 
for the auto-correlation length in traditional CT-HYB is 
given by the quotient of the number of operator pairs 
(/c/2) over the acceptance rate for removal in partition 


function space r lem z 


OH 


Nf.orr.Z 


(fc/2) 


^rem,Z 


(26) 


In principle, a similar estimate holds for the Green 
function sampling C G (n ). However, another possibility to 
arrive at an uncorrelated worm is to remove one worm 
and insert a new worm into the local trace at another 
location. If the acceptance rate for removal of a worm 
pair is r re m,W; this gives another estimate for the auto¬ 
correlation length in worm space: 

N CO rr,W ~ -, (27) 

^rem ,W 















which we employ in practice. 

It is still necessary to modify the approximations in 
Eq. (26) and Eq. (27) by the percentage of worm steps 
proposed and the percentage of hybridization operator 
steps proposed, since our new system has two different 
types of moves instead of one. We observe that the ac¬ 
ceptance rate of worm inserts and worm removals is in 
general lower when inserting four operators at once, as is 
the case for the two-particle Green’s function G^. While 
we are able to alleviate this problem partially by adjust¬ 
ing r/ 2 ), the acceptance rate is still lower due to quantum 
number violations. The reduced acceptance rate directly 
translates to an increased auto-correlation length of the 
two-particle Green’s function. 


V. ATOMIC LIMIT RESULTS 

As a first test and validation of the worm algorithm we 
consider the atomic limit. We distinguish two scenarios 
with a divergent ratio of Coulomb repulsion to hybridiza¬ 
tion strength U/V —> oo. (i) The actual atomic limit de¬ 
fined as V —> 0 for finite U, i.e., we decouple the impurity 
from the bath. In this scenario, we are still able to choose 
U freely. This allows us to control the time evolution in 
the local trace. We observe that the Green’s function 
estimators of partition function sampling fail completely 
in this case due to the absence of the hybridization func¬ 
tion. In the second scenario (ii), we keep V fixed and in¬ 
crease the Coulomb repulsion U —> oo. While the Green’s 
function estimator of partition function sampling is still 
capable of producing results for large U due to the pres¬ 
ence of the hybridization function, we observe systematic 
deviations of the error bars around r = (3/2. 


A. Atomic limit V —> 0 

The one-particle Green’s function G^ and the two- 
particle Green’s function G^ are known analytically 
in the atomic limit. On the other hand, estimators of 
the type Eq. fail completely since the impurity is 
no longer coupled to the bath. That is, measuring the 
Green’s functions by cutting hybridization lines in CT- 
HYB is no longer possible due to the absence of the hy¬ 
bridization function. The worm algorithm, on the other 
hand, is not limited by the hybridization function, as 
operators are inserted locally. As a result, the worm al¬ 
gorithm is capable of reproducing the atomic limit. 

While sampling the atomic limit with QMC algorithms 
is mainly of academic interest, we can use the analytic 
results for benchmarking. Fig. [6] shows the Green’s func¬ 
tion in the atomic limit, i.e., for an isolated impurity, 
comparing the worm algorithm and the analytic expres¬ 
sion. 

Let us now turn our focus towards two-particle quanti¬ 
ties. The measurement of four worm operators in imag¬ 
inary time is Fourier transformed into Matsubara fre- 
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FIG. 6: One-particle Green’s function G (1 \iv) in discrete 
Matsubara frequencies iv for the atomic limit of the single¬ 
orbital AIM at inverse temperature j3 = 5/D, Coulomb re¬ 
pulsion U = 1.0D and /r = 0.5U (half-filling). The balancing 
parameter was set to r / 1 * = 0.7. In the absence of any hy¬ 
bridization function, the worm algorithm (green triangles) is 
able to reproduce the analytic result (red line). Conventional 
CT-HYB is not possible. 


quencies using the particle-hole convention. The two- 
particle Green’s function G^ (iv,iv',ioS) in the particle- 
hole convention is a function of two fermionic Matsubara 
frequencies iv, iv' and one bosonic Matsubara frequency 
iui. In order to quantify results, we analyze slices of the 
full two-particle Green’s function by setting the second 
fermionic frequency to v' = ir/(3 and the bosonic fre¬ 
quency to w = 0. For comparison, we construct the 
analytic atomic limit results of the two-particle Gr een’s 
function from the expressions of the reducible vertexP^H 
A more complete discussion of the general properties of 
two-particle quantities can be found elsewhere! 3 ' Fig. Pt] 
shows the G/^(iv, n/fi, 0) slice measured using worm 
sampling and compared to the analytic result. 

We conclude that the absence of the hybridization 
function in the atomic limit results in a complete break¬ 
down of the one- and two-particle Green’s function esti¬ 
mator in partition function sampling. In contrast worm 
sampling works very well and correctly reproduces the 
analytic result for the atomic limit. 


B. Strong interaction limit U —> oo 


In principle, CT-QMC algorithms are used for interme¬ 
diate parameter ranges, but not the atomic limit itself. 
However, the strongly insulating case with high values of 
U is of interest. While here a hybridization function is 
still present for a finite bandwidth, the local time evolu¬ 
tion suppresses most of the hopping from and onto the 
impurity. 

Fig. y shows the one-particle Green’s function G^(r) 
with error bars on a logarithmic scale. Both approaches, 
partition function and worm sampling, essentially agree 
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FIG. 7: Same as Fig. [6] but now for the two-particle Green’s 
function 7T//3, 0) slice for which a balancing param¬ 

eter r= 0.155 has been employed. 



FIG. 8: One-particle Green’s function G (1 ^(t) in imaginary 
time for the single-orbital AIM with semi-elliptic conduction 
electron density of states of half-bandwidth D, U = 5.0 D and 
/r = 2.6 D (out of half-filling). The balancing parameter was 
set to r/ 1 ^ = 1.4. While the error bars of the Green’s function 
calculated by partition function sampling vanish between r £ 
[5,45] (red curve), the error bars of the Green’s function in 
worm sampling have a comparable magnitude for all values 
of r. 


for the Green’s function. However, the error bars of par¬ 
tition function sampling vanish for intermediate r-values. 
This is clearly an artifact since the error bars should be 
comparable along the whole range of r-values, as it is the 
case in worm sampling. The origin for this shortcoming is 
that hybridization pairs for intermediate r-values are no 
longer inserted, but just measured by cutting hybridiza¬ 
tion lines between operators of two operator pairs. While 
the effect on the Green’s function itself is still small, it 
already produces wrong error bars and hence maximum 
entropy spectra. Small errors may also propagate and 
get enlarged through DMFT iterations. 

We hence conclude that the worm algorithm not only 
correctly reproduces the atomic-limit but also works 
properly for large U, including error bars. As such, the 


worm algorithm provides an improvement to the conven¬ 
tional CT-HYB algorithm in the strong coupling limit. It 
also correctly reproduces the non-interacting limit mak¬ 
ing it, in principle, numerically exact over the complete 
parameter range. 


VI. TWO-PARTICLE GREEN’S FUNCTION 


In the previous section we have discussed how Green’s 
function estimators in partition function sampling lead 
to systematic errors in the absence of a hybridization 
function. This is true for any type of hybridization func¬ 
tion. Another problem arises when dealing with spin- 
orbital diagonal hybridization functions. Such a diago¬ 
nal hybridization is exact in high-symnretry cases and is 
a widely employed approximation in other systems, be¬ 
cause it mitigates the sign problem and allows for speed- 
ups due to the block-diagonalization of matrices!—*— * The 
CT-HYB algorithm is then only inserting operator pairs 
within the hybridization expansion where creation and 
annihilation operators have the same spin-orbit flavor. In 
conventional CT-HYB partition function space sampling, 
Green’s function estimator are measured by removing 
these hybridization lines. That is, one can only measure 
Green’s functions, which can be built from hybridization 
pairs with the same spin-orbit flavor. While the one- 
particle Green’s function in general fulfills this criteria 
and can be measured with such estimators (note that 
the flavor-off-diagonal one-particle Green’s function van¬ 
ishes for flavor-diagonal hybridization), this is not true 
for all components of the two-particle Green’s function. 

Especially the spin flip and pair hopping terms of the 
two-particle Green’s function are not accessible in this 
way. This is another systematic weakness of conventional 
CT-HYB partition function sampling. The worm sam¬ 
pling algorithm, on the other hand, does not suffer from 
this shortcoming. This is because four arbitrary opera¬ 
tors can be inserted into the trace. Their spin-orbit flavor 
can be chosen freely without the need to connect these 
via the hybridization function. 

In order to analyze the spin flip and pair hopping terms 
of worm sampling, we again look at the atomic limit. 
We choose the two-orbital AIM with semi-elliptic con¬ 
duction el ectro n density of states and Slater-Kanamori 
interaction ! 40 * 41 ! This local interaction includes an intra¬ 
orbital repulsion U, SU(2)-symmetric Hund’s exchange 
and pair hopping terms J, and inter-orbital interaction 
U' = U — 2 J. Fig. [9] and Fig. 10 show the spin flip term 
and the pair hopping susceptibility in the atomic limit. 
Again, we observe that worm sampling is able to repro¬ 
duce the analytic expression. 

So far we have only presented results for the spin 
flip and pair hopping term using worm sampling in the 
absence of a hybridization. While this atomic limit 
is very useful for benchmarking purposes, we are ulti¬ 
mately interested in intermediate parameters, where CT- 
QMC algorithms are predominantly used, especially for 
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FIG. 9: Two-particle spin flip Green’s function 

^iti 4 . 24 . 2 t(® z/ ’ 7r /'®’ 0) vs. iv for the atomic limit of the 
two-orbital AIM at /3 = 5/D, U = 1.0D, J = 0.4D, 
U' = 0.2D, and p = 0.5D (half-filling). The balancing 
parameter was set to r/ 2 ^ = 0.09. In the absence of any 
hybridization function, the worm algorithm (green triangles) 
is able to reproduce the analytic result (red line). 
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calculating multi-orbital systems. In order to further 
verify our results, we exploit the SU(2) symmetry of 
Slater-Kanamori-like interaction, where (S z (t)S z (0)) = 
(Sx(t)S,(0 )) holds. 

Using partition function sampling, we can calculate 
the spin susceptibility in z-direction in a straight-forward 
manner. Note that we can express S z (t) = nl(r) — nl (r) 
in terms of density operators so that (S z (t)S z ( 0)) can 
eventually be sampled by removing diagonal hybridiza¬ 
tion functions in partition function sampling. 

This is not possible for {S x (t)S x ( 0)) which is ex¬ 
pressed in terms of spin flip two-particle Green’s func¬ 
tions. While this cannot be calculated in conventional 
partition function sampling, we can do so by using 
worm sampling. Instead of looking at the imaginary- 
time resolved spin susceptibility, we verify the SU(2)- 
symmetry for the local spin susceptibility in terms of its 
Fourier transform to Matsubara frequencies Xi oc (*w) = 



0) 


FIG. 11: Local spin susceptibility ReXioc(*w) of the two- 
orbital AIM as a function of the bosonic Matsubara fre¬ 
quency iuj. Parameters: identical semi-elliptic bands of half¬ 
bandwidth D, P = 5/D, U = 1.0D, J = 0.4D, U' = 0.2D, 
and p = 0.5D (half-filling). The balancing parameter was 
set to r= 0.08. The SU(2) symmetry is conserved, as the 
S X S X susceptibility of the worm algorithm (green error bars) 
agrees well with the S Z S Z susceptibility of partition function 
sampling (red line) and worm algorithm (blue error bars). 


dT e -^(s z(x) (T)s z(x) m. 

Fig. |n| shows the spin-susceptibilities for the two- 
orbital AIM on a Bethe lattice. The worm sampling es¬ 
timate for the S X S X susceptibility in x-direction agrees 
with the S Z S Z susceptibility in z-direction, which can 
be calculated both by worm and partition function sam¬ 
pling. This further demonstrates the power of worm sam¬ 
pling to calculate general Green’s functions and suscep¬ 
tibilities. 


VII. CONCLUSION 


In this work we have demonstrated how worm sampling 
provides a solution to some systematic failures of con¬ 
ventional CT-HYB algorithms. By inserting operators 
explicitly into the local trace, we decouple the Green’s 
function measurement from the hybridization function. 
This allows us to measure the one-particle and the two- 
particle Green’s functions in situations, where the hy¬ 
bridization function is vanishing. Further, we are able 
to generate off-diagonal components of the two-particle 
Green’s function (spin flip and pair hopping terms). We 
have verified the algorithm by testing the atomic limit 
and showing the SU(2) symmetry for a two-orbital Bethe 
model. The worm algorithm supplements the hybridiza¬ 
tion expansion CT-QMC solver with a numerically exact 
procedure for estimating two-particle correlation func¬ 
tions. 
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